/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/

#include "linal/common.hxx"
#include "linal/util.hxx"

#include "linal/format.hxx"

#include "linal/operation.hxx"

#ifdef LINAL_DISPATCHER_ENABLE
#include "linal/device/operation.hxx"
#include "linal/device/device.hxx"
#include "linal/device/cache.hxx"
#include "linal/dispatcher.hxx"
#endif

#include "linal/device/host/operation.hxx"



namespace linal {
  static int ldl_nopiv_l(Hier_ A);
  static int ldl_nopiv_u(Hier_ A);

  /*!
    LDL^H factorization.
  */
  int ldl_nopiv_l(Hier_ A) {
    FLA_Obj ATL,   ATR,      A00,  A01,  A02,
            ABL,   ABR,      A10,  A11,  A12,
                             A20,  A21,  A22;

    FLA_Part_2x2( ~A, &ATL, &ATR,
                      &ABL, &ABR, 0, 0, FLA_TL );

    while ( FLA_Obj_length( ATL ) < FLA_Obj_length( ~A ) ){
      FLA_Repart_2x2_to_3x3( ATL, ATR,     &A00, &A01, &A02,
                                           &A10, &A11, &A12,
                             ABL, ABR,     &A20, &A21, &A22,
                             1, 1, FLA_BR );
      //------------------------------------------------------------
      Hier_ 
        A_11(A11), A_12(A12),
        A_21(A21), A_22(A22);

      Flat_ 
        F_11(((FLA_Obj*)A_11.get_buffer())[0]),
        F_12(~(A_12.flat())), F_21(~(A_21.flat()));

#ifdef LINAL_DISPATCHER_ENABLE
      get_dispatcher()->ldl_nopiv( FLA_LOWER_TRIANGULAR, ~F_11 );
#else
      host::internal::ldl_nopiv( FLA_LOWER_TRIANGULAR, ~F_11 );

      // ** upper(A11) := D_11 L_TL^T
      FLA_Hermitianize( FLA_LOWER_TRIANGULAR, ~F_11 );
      host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE, 
                                       FLA_UPPER_TRIANGULAR, ~F_11 );

#endif

      // ** update 
      if ( F_12.get_m() && F_12.get_n() ) {

        LINAL_LDL_L_TRSM(A_11, A_21);
        FLA_Copyt( FLA_CONJ_TRANSPOSE, ~F_21, ~F_12 );              
        host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE, ~F_11, ~F_12 );
        LINAL_LDL_L_GEMM(A_21, A_12, A_22);

      }
      
      //------------------------------------------------------------
      FLA_Cont_with_3x3_to_2x2( &ATL,  &ATR,  A00, A01, A02,
                                              A10, A11, A12,
                                &ABL, &ABR,   A20, A21, A22,
                                FLA_TL );
    }

    return true;
  }
  int ldl_nopiv_u(Hier_ A) {
    LINAL_ERROR(false, LINAL_ERROR_NOT_YET_IMPLEMENTED);
    return true;
  }
    
  int ldl_nopiv(int uplo, Hier_ A) {
    if ( !A.get_m() || !A.get_n() ) return true;

    LINAL_ERROR_DETAIL( check_uplo     ( uplo  ),
                        "%s\n %s %d\n",
                        LINAL_ERROR_CHECK_UPLO,
                        "uplo = ", uplo );

    LINAL_ERROR_DETAIL( A.get_m() == A.get_n(),
                        "%s\n %s %d\n %s %d\n",
                        LINAL_ERROR_A_SQUARE_MATRIX,
                        "length = ", A.get_m(),
                        "width  = ", A.get_n() );

    switch ( uplo ) {
    case FLA_UPPER_TRIANGULAR:        
      ldl_nopiv_u( A );
      break;
    case FLA_LOWER_TRIANGULAR:        
      ldl_nopiv_l( A );
      break;
    }

    return true;
  }
}
